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I .  INTRODUCTION 


I  - 

Two  important  parameters  of  interest  in  the  field  of  aero¬ 
dynamics  of  airfoils  are  lift  and  drag.  These  can  be  evalu¬ 
ated  either  experimentally  or  theoretically.  The  desire  for 
computational  methods  to  aid  the  design  process  is  promoted 
by  the  need  to  reduce  the  number  and  cost  of  wind  tunnel  tests. 

This  thesis  treats  the  problem  of  incompressible,  two- 
dimensional  steady  flow  about  airfoils 'or  airfoil  combinations 
at  large  angles  of  attack.  Such  flows  exhibit  strong  viscous 
flow  effects  including  regions  of  flow  separation.  Therefore 
methods  are  required  which  can  account  for  these  effects. 

Currently  there  exist  two  main  methods,  namely  the  direct 
computation  of  viscous  flows  by  means  of  the  Navier-Stokes 
equations  or  the  so-called  viscous/inviscid  interaction  method. 
The  former  approach  is  more  straightforward  but  also  much  more 
expensive  and  time-consuming.  Therefore,  the  latter  approach 
is  to  be  preferred  if  it  can  be  shown  that  it  produces  good 
agreement  with  the  available  experimental  results. 

It  is  the  purpose  of  this  thesis  to  contribute  to  the 
evaluation  of  the  viscous/inviscid  interaction  method.  To 
this  end,  the  viscous/inviscid  computer  codes  developed  by 
Cebeci  and  collaborators  at  the  Douglas  Aircraft  Company  were 
obtained  and  applied  to  several  airfoils. 

In  addition,  a  separate  panel  method  was  formulated  and  pro 
grammed  in  order  to  obtain  the  inviscid  flow  over  airfoil 


combinations . 


The  basic  equations  are  formulated  in  Chapter  II.  Chapter 
III  addresses  the  problem  of  inviscid  flow  calculations  using 
the  panel  method.  In  Chapter  IV  the  solution  of  the  boundary 
layer  equations  by  means  of  the  Cebeci-Keller  box  method  is 
explained.  Finally,  Chapter  V  describes  the  viscous/inviscid 
interaction  problem  and  presents  results  of  computations. 


II 


A.  INTRODUCTION 

In  this  chapter,  the 


exact  so. 
classif  i< 


r.liiptic 


/  «■  0  "  u  *  O  * 


Meaning 


!  Example 
!  Equations 


Upstream 

Influence 

No  Upstream 
Influence 

•  Laplace 

•  Steady  Navier- 
Stokes 

•  Thin  Shear 
Layer 

For  example,  the  Laplace  equation  is  elliptic.  Its  solution 
would  have  to  be  repeated  for  many  iterations  so  that  the  up¬ 
stream  influence  can  be  gradually  propagated  (panel  method  in 
Chapter  III) ,  but  the  thin  shear  layer  equations  are  parabolic. 
Their  numerical  solution  can  be  obtained  by  marching  in  the 
downstream  direction  only  (Box  method  in  Chapter  IV)  . 


B.  DERIVATION  OF  GENERAL  EQUATIONS 

The  continuity  equation  and  the  Navier-Stokes  equations 
are  basic  for  an  aerodynamic  analysis.  We  start  with  the  basic 
physical  concepts  and  derive  the  general  equations  for  2-D, 
unsteady,  compressible,  viscous  flow. 

1 .  Continuity  Equation 

One  of  the  basic  laws  of  "Newtonian  mechanics"  states 
that  mass  can  neither  be  created  nor  destroyed.  Therefore, 
for  a  fixed  control  volume  (see  Figure  2.1),  the  principle  of 
mass  conservation  can  be  stated  that  the  net  mass  flow  rate 
into  and  out  of  the  control  volume  equals  the  time  rate  of 
change  of  mass  within  the  control  volume. 

If  the  central  point  'P'  has  representative  fluid  proper 
ties  (velocity,  density,  pressure,  etc.),  then  properties  at 
other  locations  can  be  obtained  by  Taylor  series  expansions. 
Therefore  the  x-component  of  the  velocity  at  the  center  of  the 
positive  x-face  (right-hand  face)  is 


Figure  2.1.  Control  Volume  for  2-D 


As  dx  goes  to  zero,  all  higher  order  terms  will  be  dropped, 
so  that  only  the  first  two  terms  will  be  considered.  Similarly, 
the  density  at  the  center  of  the  positive  x-face  is 


P 


+ 


9  p  dx 
9x  2 


(2.2) 


Then  the  mass  flow  rate  out  of  the  positive  x-face  is 


Mass  Flow  Rate  Out 


(Density) (Velocity) (Area) 


••><»♦&¥ ♦ 


[OU  +  !j(pu)3*ldy 


. . . ) dy 

(2.3) 


By  the  same  procedure  the  mass  flow  rate  into  the  control  volume 
through  the  negative  x-face  (left-hand  side)  is 


Mass  Flow  Rate  In 


(2.4) 


=  [pu  -  |j(pu)  ^]dy 


From  Eqs.  (2.3)  and  (2.4),  we  get  the  net  mass  flow  rate 
through  the  control  volume  in  the  x-direction. 


Net  Mass  Flow  Rate  = 


[pu  -  |^(pu)^] dy  -  [pu  +  |^(pu)^r]dy 


=  -  ^-(pu)dx  dy 


(2.5) 


In  a  similar  manner,  the  net  mass  flow  rate  in  the  y-direction 


^y(pv)dx  dy 


(2.6) 


The  total  mass  flow  rate  through  the  control  volume  is 
obtained  by  summing  Eqs.  (2.5)  and  (2.6). 


Total  Net  Mass  Flow  Rate  = 


-l|^(pu)  +|^-(pv)  ]  dx  dy  (2.7) 


Next,  we  consider  the  time  rate  of  change  of  mass  within 
the  control  volume. 


Time  Rate  of  Change 
of  Mass 


=  7-r-(p  dxdy) 


Now  we  combine  Eqs.  (2.7)  and  (2.8)  using  the  concept  of 
conservation  of  mass.  Then 


tfx(pu)  +  f^pv) idx  dy 


|£dxdy 


(2.9) 


Therefore  we  obtain  the  general  form  of  the  continuity 
equation  for  two-dimensional  flow  as 


9 (pu)  9 (pv)  ,  9  P  _  n 

ax  9y  at 


(2.10) 


For  steady  or  unsteady  incompressible  flow,  Eq.  (2.10)  reduces  to 


3u  3v 
3x  ay 


(2.10a) 


2 .  Navier-Stokes  Equations 

Newton's  second  law  of  motion,  when  mass  is  conserved, 
equivalently  states  that  the  rate  of  change  of  the  momentum  of 
a  body  equals  the  sum  of  the  forces  applied  to  that  body,  or 


l  F  =  hF^V) 


(2.11) 


In  considering  a  small  volume  element  of  fluid,  there  are  two 
types  of  forces  to  be  considered,  namely  surface  forces  which 
are  acting  on  the  surface  and  body  forces  which  are  acting  on 
the  fluid  inside  the  elemental  volume,  such  as  gravity  (see 
Figure  2.2). 


*  rj?r; 


y 

* 


-■>  x 


1  3avx  dx 

■  ayx  +  ~3x  2 


Gravi  ty 


3oxx  dx 
axx  +  3x  2 


V 


3a  , 

a - Z* 

yx  3y  2 


3a  , 

a _ 

yy  »y  2 


Figure  2.2.  Forces  Acting  on  the  Fluid 


Assuming  that  the  stresses  are  known  at  point  1 P’,  we  get 
expressions  for  the  stresses  on  the  fluid  element  surfaces 
by  Taylor  series  expansion. 


Net  force  in 
x-direction  due 
to  Normal  Stresses 


3a 

r  .  XX 

[a  +— =r — 

1  xx  3x 


3a 


xx 

dX 


3a 


XX 


3x 


dx  dy 


(2.12) 


Net  force  in 
x-direction  due 
to  Shear  Stresses 


3a  ,  3a  . 

|o  +_JI2&dx-[a  -JE'Sijdx 

1  yx  dy  2  1  yx  oy  2 


3a. 


3y 


^  dx  dy 


(2.13) 


Therefore  the  total  surface  forces  in  the  x-direction  are 
formed  by  summation  of  Eos.  (2.12)  and  (2.13). 
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I  fx  (surface)  =  l-jf2  +  -J2S]  dxdy 


(2.14) 


Let  f  (body)  be  defined  as  the  body  force  per  unit  mass  with 
the  following  components : 


t  (body)  =  X  i  +  Y  j 


Thus,  the  x-component  of  the  body  force  is 


f  (body)  =  m  X 

X 


=  p  dxdy-1-x 


(2.15) 


Adding  Eqs.  (2.14)  and  (2.15)  provides  the  total  force  in  the 
x-direction . 


f  (surface)  and  f  (body) 

X  X 


9a  9a 

+  T?  +  Tf  1  «xdy 


(2.16) 


Now  we  consider  the  rate  of  change  of  the  momentum  of  the  fluid. 

Let  us  take  the  x-component  only.  Then,  since  the  mass  is  constant. 


^(mV) 


.  ,  ,  9u  ,  9u  .  9u, 
p  dxdy  (u  —  +  v  ^  ^ 


(2.17; 


*-'v  • 


because,  u  =  f [x (t) ,y (t) , t] ,  and  by  chain-rule 


du  _  3u  dx  3u  dy  3u  3t 
dt  3x  dt  3y  dt  3t  3t 


u 


3u 

3x 


+  v 


3u  3u 
3y  9t 


Substitution  of  Eqs.  (2.16)  and  (2.17)  into  the  x-component 
of  Eq.  (2.11)  gives  the  final  equation  of  motion. 


pX  + 


3a 


xx 


3x 


r  3u 
0[u  35? 


+  V 


(2.18) 


Now  we  want  to  show  the  stress  in  terms  of  the  velocity  com¬ 
ponents.  In  this  thesis  we  will  consider  only  simple  "Newtonian 
fluids  obeying  Stokes'  law.  This  means  that  the  'extra'  stress 
(above  hydrostatic  pressure)  is  proportional  to  the  rate  of 
strain.  With  the  definition. 


Extra  Stress  =  constant  x (rate  of  strain) 


and  introducing  p  =  coefficient  of  viscosity, 


xx 


-p  +  ^  <&> 


xy 


,3u  ,  3v, 
u(3?  +  35?’ 


where  the  pressure  in  an  incompressible  fluid  is  seen  to  be 
equal  to  minus  one-third  the  sum  of  the  three  normal-stress 


2] 


components  in  view  of  Eq.  (2.10a).  In  two  dimensions  then, 


o  +  a  =  -2P  .  Eq.  (2.18)  then  becomes,  if  the  body  force 
xx  yy 

is  neglected 


3u 

3t 


+ 


u 


V 


3u 

3y 


1  3P  1  9axx  1  9axy 

p  3x  p  3x  p  3y 


(2.19) 


Substitution  then  produces,  for  incompressible  flow. 


u 


v 


3u 

3y 


3x 


3y 


(2.20) 


where  v  =  p/p  =  kinematic  viscosity,  and  similarly  for  the 
y-component 


3  v 
3 1 


+ 


3v 

3 


v 


3v 

3y 


-  i  |§  + 

p3y  ~  2  ~  2 


3x 


3y 


(2.21) 


These  are  the  well-known  Navier-Stokes  equations  for  two- 
dimensional  incompressible  viscous  flow. 


C.  INVISCID  FLOW  EQUATION 

All  real  fluid  flows  are  viscous,  but  inviscid  flow  can 
be  assumed  outside  of  a  thin  boundary  layer  and  a  narrow  wake 
behind  the  body  for  large  Reynolds  numbers.  This  is  the  reason 
why  the  inviscid  flow  equations  are  important  even  though  they 
represent  an  ideal  case.  If  the  flow  far  upstream  is  uniform 
then  it  is  also  irrotational .  This  allows  the  introduction  of 


the  velocity  potential. 


1 .  Velocity  Potential 

The  velocity  potential  <j>  is  a  function  whose  gradient 
represents  the  fluid  velocity.  Thus 


3<t> 

3x 


u 


(2.22) 


where 


0  =  <J>  (x,y,  t) 

Therefore,  the  significance  of  the  velocity  potential  lies 
in  the  fact  that  one  equation  for  <J>  can  be  used  rather  than 
three  equations  for  the  velocity  components. 

2 .  Laplace  Equation 

For  steady,  incompressible  flow,  the  continuity  equa¬ 
tion  (2.10)  becomes 


+  £  ■  0 

This  equation  can  be  written  in  terms  of  velocity 
potential  <j>  by  substituting  Eq.  (2.22).  Thus 

^4  +  i4  =  0  (2.24) 

3x  3y 

2 

This  is  the  well-known  Laplace  equation  (vector  form  is  V  4>  =  0)  . 
It  is  a  linear  equation  which  makes  it  possible  to  apply  the 
principle  of  linear  superposition.  For  instance. 
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If  4>^,  <J>2,  .  ..,  4>n  satisfy  V  <p  =  0,  then  also 

2 

<J>  =  +  <P0  +  ...  +  <j>  satisfies  V  <p  =  0 . 


Thus  the  flow  resulting  from  the  superposition  of  incompressi¬ 
ble,  irrotational  flows  is  also  an  incompressible  and  irro- 
tational  flow.  This  superposition  principle  makes  it  possible 
to  build  up  quite  complicated  flow  from  a  few  simple  solutions 
of  Laplace's  equation.  The  singularity  (or  panel)  methods 
presented  in  the  next  chapter  are  based  on  this  idea. 


D.  THIN  SHEAR  LAYER  EQUATIONS 

High  Reynolds  number  flows  over  airfoils  (and  other  con¬ 
figurations)  generate  a  very  thin  shear  layer  (boundary  layer) 
close  to  the  airfoil  surface.  If  6  denotes  the  boundary  layer 
thickness  and  x  the  distance  from  the  leading  edge  of  a  flat 
plate,  then  it  is  well-known  that 


<5  (x)  ~  /vx/U 


6  (x) 


-  /1/Re 


where  v  is  the  kinematic  viscosity  and  Re^  =  Ux/v.  This  formula 
shows  that 


5  (x) 


<<  1 


>>  1 


Hence  the  flow  outside  of  the  boundary  layer  can  be  considered 
to  be  inviscid,  but  the  effect  of  viscosity  cannot  be  neglected 
within  this  layer.  Nevertheless,  the  Navier-Stokes  equations 
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for  steady  incompressible  flow  can  be  simplified  because  of 
the  fact  that  6  is  much  smaller  than  the  representative  length 
of  the  airfoil  (the  chord  length) .  This  can  be  deduced  from 
the  Navier-Stokes  equations  as  follows: 


1 

P 


3P 

3x 


v 


3u 

3y 


(2.25) 


1  3P 
P  3y 


+ 


u 


v 


3v 

3y 


(2.26) 


u  is  now  replaced  by  a 
x  is  now  replaced  by  a 
y  is  now  replaced  by  a 

3  u 

Then  can  be  expressed 
^  can  be  expressed 

^  can  be  expressed 


typical  value,  say  Ue ; 
typical  value,  say  l; 
typical  value,  say  6. 

by  Ue/<5; 
by  Ue/!l; 

by  pu2/& 
e 


(because  P  and  Ue  are  related  by  the  Bernoulli  equations) . 
Therefore  the  x-component  terms  of  the  Navier-Stokes  equation 
can  be  estimated  to  have  the  following  magnitudes: 


The  magnitude  of  the  term  v 

of  the  continuity  equation  + 
5U  dx 

v  -  and  therefore 


follows  from  the  application 


i .e . , 


0  or 


v 


3u 

3y 


6 


l  8 


The  two  viscous  terms  are  of  vastly  different  magnitude 

because  32u/3x2  -  Ue/£2  and  32u/3y2  -  U^/6  2  hence 

2  2  '  2  2  2  2 
3  u/3y  >>  3  u/3x  and  3  u/3x  can  therefore  be  neglected 

2  2  2  2 
compared  to  3  u/3y  .  Finally,  the  term  v  3  u/3y  must  be  of 

the  same  magnitude  as  the  other  terms  if  the  influence  of 

viscosity  is  to  be  retained.  The  y-component  terms  of  the 

Navier-Stokes  equations  are  easily  estimated  to  be  smaller 

than  the  x-component  terms  because 


u 


3v 

3x 


U 


_6_ 

,2 


6 

and  hence  are  smaller  by  a  factor  j-.  Therefore  the  two 
equations  reduce  to 


u 


3u  ,  3u 
3?  +  v  37 


1  3P 

p  3x 


+  v 


(2.27) 


P  =  0  (2.28) 

3y 


By  adding  the  continuity  equation  (Eq.  (2.23))  to  these  rela¬ 
tions,  we  get  the  basic  equations  to  describe  laminar  flow 
thin  shear  layers. 


E .  TURBULENT  FLOW 


We  must  deal  with  the  instantaneous  properties  in  the 
turbulent  flow.  Thus  the  time-mean  value  concept  is  applied: 

u  =  u  +  u ' 

where  u  is  the  time-mean  value,  and  u'  is  the  fluctuating 
portion.  Similarly, 


v  =  v  +  v 1 


P  =  P  +  P' 


Introducing  these  relations  into  Eq.  (2.20) 


1  3P 
p  9x 


+  v( 


9u 1 u 1  _  9u 1 v 1 
3x  9y 


(2.29) 


We  can  see  that  pu'u'  and  pu ' v '  correspond  to  a  normal  stress 
and  a  shear  stress.  We  call  these  terms  turbulent  stresses 
or  Reynolds  stresses. 

Similar  analyses  can  be  done  for  the  y-component  equations 
and  z-component  equations  in  the  three-dimensional  case.  The 
extra  Reynolds  stresses  can  be  summarized  by  the  following 


array. 


III.  PANEL  METHOD 


A.  INTRODUCTION 

The  panel  method  was  developed  in  the  1960's  at  McDonnell 
Douglas  by  Smith  and  Hess  as  a  numerical  approach  for  2-D 
and  3-D  potential  flow  problems.  This  chapter  presents  the 
application  of  the  panel  method  to  2-D  problems  about  one  or 
two  bodies.  The  basic  idea  consists  in  representing  the  flow 
past  a  body  by  a  distribution  of  singularities  (sources, 
sinks,  vortices)  on  the  surface  such  that  the  body  surface 
becomes  a  streamline. 

The  numerical  approach  requires  some  approximations  (the 
assumption  given  in  parentheses  refers  to  our  approach) . 

A.  The  surface  of  the  body  is  replaced  by  a  finite 
number  of  elements  (straight-line-elements) . 

B.  The  condition  of  tangential  flow  is  satisfied  at  a 
finite  number  of  points,  the  so-called  control- 
points  (midpoints  of  elements) . 

C.  The  singularity  distribution  of  each  element  is 
approximated  by  some  kind  of  analytical  functions 
(singularity  strengths  are  assumed  to  be  constant 
along  any  one  element,  but  vary  from  element  to 
element) . 

The  advantages  of  the  panel  method  in  comparison  with  other 
procedures  are: 

A.  The  panel  method  does  not  include  an  approximation  in 
the  physics — thin  airfoil  theory  does. 

B.  The  panel  method  can  be  easily  applied  to  both  2-D 
and  3-D  problems — a  virtually  unsolvable  task  for 
conformal  mapping  procedures,  which  are  confined 
to  2-D  configurations  only. 
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C.  The  panel  technique  can  be  readily  extended  to  flow 
fields  past  several  bodies — a  task  which  causes  at 
least  some  troubles  in  "classical"  mapping  techniques. 

The  method's  versatility  has  been  proven  in  various  extensions, 
e.g.,  hydrofoils,  cascades,  nozzles,  and  even  complete  air¬ 
craft*  Since  its  origin  the  method  has  been  improved  by 

features  like  higher  order  approximations  to  both  body  surface 
and  singularity  distribution,  taking  account  for  compressi¬ 
bility  effects,  and  inclusion  of  wake  models.  Today  the  panel 
method  is  probably  our  most  powerful  tool  in  analyzing  poten¬ 
tial  flows. 

B.  NONLIFTING  FLOW  PAST  A  BODY 

The  effects  of  lift  (respectively,  camber  and  angle  of 
attack)  and  displacement  (resp.  thickness)  can  be  studied 
separately  because  of  the  linearity  of  Laplace's  equation. 

This  section  is  concerned  about  displacement  flows  due  to  the 
thickness  of  bodies,  a  flow  which  is  usually  represented  by 
sources  and  sinks. 

We  will  first  draw  the  reader's  attention  to  a  single 
straight-line-element,  along  which  sources  of  constant  strength 
are  distributed.  This  simple  case  allows  us  to  explain  the 
basics  of  the  panel  technique.  The  source  strength  A  is  de¬ 
fined  as  the  volume  of  fluid  discharged  per  unit  area.  Since 
fluid  is  ejected  perpendicular  to  the  panel's  surface  in  both 
directions,  discharge  velocities  are  half  of  the  source 
strength.  The  boundary  condition  for  an  inclined  panel  requires 


that  the  normal  component  of  the  free  stream  velocity  is 
balanced  by  this  discharge  velocity  (see  Figure  3.1). 


X 

2 


-  V  cos  6 
00 


(3.1) 


panel 


Figure  3.1.  Boundary  Condition  at  One  Inclined  Panel 

This  relation  between  a  geometric  quantity  (£)  and  the  unknown 
source  strength  establishes  tangential  flow  on  the  panel 
surface . 

Things,  which  are  obvious  for  a  single  panel,  become 
slightly  more  complicated  for  a  structure  of  panels.  Mutual 
interference  of  source  panels  i.e.,  each  panel  induces  a 
velocity  at  other  panels,  causes  the  complication.  While  the 
boundary  condition  of  a  single  panel  had  been  set  up  by  glancing 
at  a  simple  geometry  sketch,  we  now  better  switch  to  a  syste¬ 
matic  procedure,  emphasizing  the  concepts  of  velocity  poten¬ 
tial  and  superposition.  We  consider  a  2-D  closed  body, 


approximated  by  several  panels  and  inclined  at  an  angle  to 
the  oncoming  flow.  Our  goal  is  to  derive  a  relation  for  the 
unknown  source  strengths  from  the  condition  of  tangential  flow 
at  the  control  points.  To  this  end  we  will  give  the  velocity 
potentials  of  a  single  source,  a  source  panel,  and  a  closed 
body  built  up  by  a  source  panel,  in  that  order. 

Radial  streamlines  and  concentric  circular  equipotentials 
characterize  one  of  the  very  basic  potential  flows,  the  single 
source  flow.  Its  velocity  potential  is  defined  by 

single  source  ( x , y )  =  ^  Unr  (3.2) 

where  r  is  the  distance  from  (x,y)  to  the  source.  Arranging 
single  sources  on  a  straight  line  element  corresponds  in 
terms  of  the  velocity  potential  to  a  summation  of  single  poten¬ 
tials,  in  the  limit  of  an  infinite  number  of  sources  to  an 
integration  over  the  panel  length.  Thus  the  velocity  potential 
of  a  source  panel  can  be  written  as 

A  £ 

source  panel  0(x,y)  =  j—  j  in  r  dS  (3.3) 

m  source  panels,  representing  a  body,  induce  a  flow  field, 
whose  velocity  potential  at  any  point  (x,y)  is  given  by 

,  l  . 

m  m  A  .  j 

=  l  4>-j  (x,y)  =  l  f  in  r  dS . 


0 (x,y) 


(3.4) 


Figure  3.2.  Designations  for  Calculation 

We  call  <J>  the  potential  of  the  flow  disturbances  due  to  the  di 
placement  flow.  The  total  velocity  potential  of  a  nonlifting 
flow  results  from  a  superposition  of  this  displacement  flow 
and  a  uniform  flow,  which  is  inclined  at  an  angle  a  to  the 
x-axis . 


,  I . 
m  A  .  j 

$>(x,y)  =  v  cosax  +  V  sin  ay  +  £  —1  f  £n  r .  dS  . 

j=l  ^  0  D  J 


Recalling  the  definition  of  the  velocity  potential  (velocity 
equals  gradient  of  potential) ,  the  boundary  condition  of  tan¬ 
gential  flow  takes  the  form 


3  & 

-r —  =0  on  the  surface 

on 


Applying  this  condition  within  the  framework  of  the  panel 
method  provides  a  system  of  equations, 
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,  l . 
m  A.  d 


V  cosa  cos 8 •  +  V  sina  sing  +  V  _1  J  __(£n  r.  .)ds. 
00  1  00  1  ztt  '  n  dn .  11  1 

1=1  0  i  J 


for  i  =  1 , . . .  ,m 


(3.5) 


which  establishes  zero  normal  velocity  at  all  control  points. 
This  linear  system  can  be  solved  for  the  unknown  source 
strengths  after  the  integrals  have  been  evaluated. 

Concept  of  influence  coefficients 

The  influence  coefficient,  defined  as  the  normal 

velocity  at  the  i  panel  due  to  a  source  distribution  of 
2Tr-strength  at  the  j*"*1  panel. 


/.  3n“^n  rij>  dSj 
j  1  J  J 


(3.6) 


The  contribution  to  the  normal  velocity  at  the  ith  panel  by 

i.  u 

the  actual  source  distribution  of  the  j  panel  is  the  product 
of  Aj/2fr  and  the  influence  coefficient  I^j  .  To  compute  the 
influence  coefficient  we  must  substitute 


r 


ij 


+ 


2 


in  Eq.  (3.6)  and  carry  out  the  differentiation. 
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where 
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we  obtain 


1  j  lxM  ~  +S.cos0.)]cos6i+  [yM  -  (yB  +S_.  sine  )  ]  sinBj^ 

.  =  j  - i - 3 _ _i _ 3 _ _  _ 

lj  0  [  xM  -(xB  +Sjcos0j)]2+[yM  -(yB  +SjSin0  J]2 


(3.8) 


Equation  (3.8)  is  expressed  in  terms  of  only  and,  after 


some  manipulations,  the  integral  may  be  written  in  the  form 
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where 


5  =  (V'V)COS  +  (^M."yB.,Sin 


i  3 


l  3 


c  =  cos  0 . cos  6.  +  sin  0 . sin  B- 
3i  3i 
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e 


=  2  cos0,(xM  -xB  )  +  2  sine  (y  -yB  ) 

J  i  j  J  i  j 


f  =  (XM._XB.)2  +  {yM.-yB.)2 

i  d  id 


The  integration  can  now  be  performed  analytically. 
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(3.10) 


Determination  of  unknown  source  strengths  A 

Equation  (3.5),  expressed  in  terms  of  Eq.  (3.9)  and  divided 
by  V  ,  takes  the  form 
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(3.11) 


where 
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or  in  the  more  convenient  matrix  form 
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The  above  set  of  linear  equations  can  be  solved  for  Aj  by 
Gauss'  elimination  method  or  any  other  linear  equation 
algorithm. 

On-body  velocities 

i.U 

The  velocity  at  the  midpoint  of  the  i  panel,  Vg  can 

1 

be  obtained  by  a  spatial  derivative  of  the  velocity  potential 
in  tangential  direction. 


3<P  (xi,yi) 


3x .  3y .  m  3  (Jin  r .  .  ) 
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where 
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(3.14) 


where  J ^.  =  /  (Jin  r^.)dS.  denotes  the  tangential  velocity 

1-)  j  l  13  3 

at  the  ifc^  panel  due  to  a  source  distribution  of  27T-strength 
at  the  j panel. 


The  calculation  of 
of  Iij,  so 


follows  the  same  procedure  as  that 
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13  13 

Positive  signs  of  on-body  velocities  indicate  that  velocities 
are  oriented  in  the  direction  of  the  surface  coordinates, 
while  negative  signs  indicate  opposite  directions  of  veloci¬ 
ties  and  surface  coordinates.  The  positive  direction  of 
surface  coordinates  is  defined  clockwise.  Therefore  positive 
values  of  on-body  velocities  are  to  be  expected  on  the  upper 
surface,  negative  values  on  the  lower  surface. 

Off-body  velocities 


Streamlines  can  be  determined  by  computing  velocities  at 
off-body  points  and  using  a  numerical  quadrature  to  progress 


from  one  point  to  another  point  on  a  streamline.  The  velocity 
components  can  be  expressed  according  to 
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Normalizing  the  above  equations  by  the  free  stream  velocity 
and  abbreviating  the  integrals  simplify  the  relations  to 
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X  V 

I^j  and  are  again  influence  coefficients,  whose  evaluation 

can  be  adopted  from  the  already  introduced  procedure. 
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Cx  =  cos  0. 
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C.  CIRCULATORY  FLOW 

While  inviscid  2-D  flow  theories  are  not  capable  of  pre¬ 
dicting  drag  characteristics,  information  about  lift  can  be 
provided  by  them.  Creation  of  lift  is  closely  related  to  a 
type  of  flow  called  circulatory  flow.  We  mentioned  already 
that  the  flow  around  a  lifting  airfoil  can  be  decomposed  into 
two  elementary  flows,  i.e.,  displacement  flow  and  circulatory 
flow.  Circulation  and  circulatory  flows  are  the  subjects  of 
this  section. 

The  early  approaches  of  airfoil  theory  emphasized  a  flow 
model  in  which  the  airfoil  was  represented  by  an  infinitely 
thin  vortex  sheet  only.  This  so-called  thin  airfoil  theory 
predicts  lift  quite  well,  because  lift  depends  primarily  on 
circulatory  flow.  Unfortunately  a  straightforward  extension 
of  vortex  sheets  to  "surface  singularity"  method  is  impossible. 
Therefore  aerodynamicists  have  proposed  a  couple  of  flow  models, 
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which  allow  the  implementation  of  circulatory  flows  in 


"surface-singularity"  methods.  Examples  are: 

(1)  Smith  and  Hess  represent  circulatory  flows  by  a  com¬ 
bination  of  source  and  vortex  distributions.  [Ref.  1] 

(2)  Martensen  prefers  vortex  distributions  only,  but  states 
the  problems  in  terms  of  the  stream  function.  [Ref.  2] 

(3)  Davenport  makes  use  of  linearly  varying  vortex 
distributions.  [Ref.  3] 

Our  approach  follows  the  ideas  of  Smith  and  Hess.  These  circu¬ 
latory  flows  are  composed  of  a  vortex  distribution,  which  is 
constant  along  all  and  for  all  panels,  and  a  source  distribution 
of  conventional  shape. 

We  start  at  the  very  beginnings  of  vortex  flows.  Concentric 
circular  streamlines  and  radial  equipotentials  characterize  the 
flow  field  of  a  single  vortex.  Its  velocity  potential  can  be 
written  as 


single  vortex  <p(x,y) 
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2tt 


arctan 


y-y. 


x-x. 


(3.20) 


with  (xv,yv)  as  the  center  of  the  vortex.  A  structure  of  m 
vortex  panels  induces  a  flow  field,  whose  velocity  potential 
at  any  point  (x,y)  is  given  by 

l  . 

m  .  ]  y-y . 

4>(x,y)  =  T  l  (-  Jj-)  I  arctan  — — 1  ds.  (3.21) 

j=l  0  j  J 

This  flow  field  differs  in  two  important  points  from  the  non¬ 
lifting  flow  field: 

(1)  It  violates  the  condition  of  tangential  flow. 


(2)  The  unknown  singularity  strength  T  cannot  be 
determined  immediately. 


The  task  of  determining  circulation  must  be  postponed  to  the 
implementation  of  the  Kutta  condition.  Temporarily  we  set 
the  vortex  strength  equal  to  one.  Tangential  flow  must  be 
established  by  the  aid  of  an  additional  source  distribution. 
Strengths  of  this  additional  source  distribution  must  be  com¬ 
puted  according  to  the  condition  that  the  normal  velocities 
due  to  the  vortex  distributions  at  the  control  points  are 
balanced  by  the  normal  velocities  due  to  the  additional  source 
distribution. 


f  | —  £n  (r  .  .  )  dS  .  = 

•  f,  2ir  J  .  3n  id  j 

J  —  J-  J  1 


m  2  .  y . -y . 

.y.  k7<tan"  jrWL>asj 

D=1  D  i  i  D  J 


Abbreviating  the  integrals  by  the  above  defined  influence 


coefficients,  we  get 


l  M 

j  =  l  11 


(!)  j(s)  = 


l  I.(v 

j=l  13 


for  i  =  1 , .  .  .  ,m 


(3.22) 


where 


X .  are  the  unknown  strengths  whose  effect  is  intended 
to  balance  the  normal  velocities  induced  by  a  unit 
vortex  distribution. 

( s ) 

I.  .  is  the  normal  influence  coefficient  due  to  a  source 
;L•,  distribution. 

lfV^  is  the  normal  influence  coefficient  due  to  a  vortex 
1-)  distribution. 

( s ) 

J.  .  is  the  tangential  influence  coefficient  due  to  a 

source  distribution. 

jfV^  is  the  tangential  influence  coefficient  due  to  a 
1-J  vortex  distribution. 


Since  influence  coefficients  of  source  and  vortex  distributions 
are  related  by  l|Y^  =  ,  the  above  equation  can  be  expressed 

according  to 


(1)  T  (s)  _ 


l  \\ ^ 

j-1  3  X3 


\  J. .  for  i  =  1 ,  . . .  ,m 


(3.23) 


By  solving  this  system  for  A  j  ,  we  determine  the  properties 

of  circulatory  flow  of  unit  strength. 

Calculation  of  disturbance  velocity  due  to  unit  circulatory 
flow 

The  disturbance  velocity,  v|v^ ,  is  composed  of  two  parts, 
one  due  to  the  constant  vortex  distribution,  the  other  due  to 
the  additional  source  distribution 


V(v) 

i 


i  j  j'v'  +  y  i!l)  j!s) 


(3.24) 


Making  use  again  of  the  relation  between  influence  coefficients 

( J fY^  =  lf®^),  we  have 
13  13 


V(v) 

1 


I  I**’  +  I  * 

i=l  13  i=i  3 


(i)  T(s) 

j  jij 


(3.25) 


D.  SYNTHESIZING  A  COMBINED  FLOW 

The  Kutta  condition  serves  as  matching  condition  for  nonlift¬ 
ing  and  circulatory  flow.  These  two  basic  flows  must  be 
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Since  the  panel  method  does  not  permit  the  evaluation  of 
velocities  at  the  trailing  edge,  the  Kutta  condition  is  satis¬ 
fied  approximately  by  requiring  that  velocities  at  the  control 
points  of  the  rearmost  panels  have  equal  magnitude.  Therefore 
the  rearmost  panels  should  be  chosen  short  so  that  flow  at 
their  midpoints  will  effectively  represent  that  at  the  trailing 
edge. 

Determination  of  circulation 

Suppose  1st  and  m*"*1  panels  are  the  closest  panels  to  the  trail 
ing  edge  on  the  lower  and  upper  surface,  then  we  can  write  the 


Kutta  condition  as 


V(N) 


+  rv 


(v) 


m 


(3.26) 


where  v\  denotes  the  tangential  velocity  in  nonlifting  flow. 

Equation  (3.26)  can  be  solved  for  the  circulation  r. 

Calculation  of  on-body  velocities  and  of  pressure  coefficients 

Three  parts  contribute  to  the  total  velocity:  free  stream, 

disturbance  due  to  displacement  flow  and  disturbance  due  to 

(N) 

lifting  flow.  Say  V  designates  the  velocity  due  to  the 
nonlifting  flow  including  the  free  stream  component  and 
represents  the  velocity  due  to  a  lifting  flow  of  unit  circu¬ 
lation.  Then  the  total  tangential  velocity  at  the  midpoint  of 
th 

the  i  panel  is  given  by 


V. 

1 


v  +  rvfv) 


(3.27) 


Once  the  velocity  has  been  computed,  the  pressure,  customarily 

expressed  by  means  of  a  dimensionless  coefficient  C  ,  is 

F 

determined  by  Bernoulli's  equation: 


(3.28) 


Addendum:  More  than  one  body  configuration. 

One  of  the  main  advantages  of  the  panel  method  is  its 
easy  extension  to  multi-element  airfoils.  As  a  matter  of 


fact  even  flow  past  an  infinite  number  of  bodies  can  be  solved 
by  means  of  the  panel  method,  if  these  bodies  are  arranged  in 
form  of  cascades.  The  minor  changes,  which  are  necessary  to 
apply  the  panel  method  to  a  finite  number  of  bodies,  include: 

(1)  The  overall  scheme  must  provide  a  circulatory  flow 
for  each  lifting  body.  (The  number  of  nonlifting 
flows  remains  one.) 

(2)  Flow  past  each  body  with  lift  is  subject  to  a  Kutta 
condition.  Accordingly  the  numbers  of  equations 
requiring  zero  load  at  the  trailing  edge  equals  the 
number  of  circulatory  flows,  which  allows  the 
definite  determination  of  each  lifting  body's 
circulation . 

Figure  3.4  illustrates  the  superposition  of  nonlifting-  and 
circulatory  flows  for  a  two  element  airfoil. 


E .  EXAMPLES 


This  section  illustrates  the  capabilities  of  the  program 
"PANEL”  which  can  be  applied  to  2-D  potential  flow  problems 
past  one  or  two  bodies. 

Flow  past  one  circular  cylinder 

The  source  panel  technique  is  applied  to  the  flow  past 
a  circular  cylinder.  This  case  is  regarded  as  nonlifting, 
i.e.,  the  cylinder  experiences  no  force  perpendicular  to  the  free 
stream.  As  sketched  in  Figure  3.5/  the  surface  of  the  cylinder 
is  approximated  by  eight  panels  of  equal  width.  For  zero  angle 
of  attack,  Eq.  (3.5)  reduces  to 

A  .  m  A.  3 { £n  r  .  .  ) 

Vg°  +  —  +  I  2 i  I.  in.  13  dsj  -  0  <3'29> 

3=1  3  i  J 

j*i 

Solving  a  set  of  8  simultaneous  algebraic  equations,  the  source 
strengths  and  the  pressure  coefficients  can  be  determined. 


*  y 


Figure  3.5.  Arrangement  of  Panels  on  a  Circular  Cylinder 


The  results  are  shown  in  Figure  3.6  where  they  are  com- 

2 

pared  with  analytical  results  (C  =  1  -4sin  6). 

r 
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Figure  3.6.  Pressure  Coefficient  on  a  Circular  Cylinder 
Obtained  by  Using  Eight  Source  Panels 
(Marked  by  0)  in  Comparison  with  the  Exact 
Solution 

This  example  demonstrates  the  power  of  the  panel  method. 
However  the  reader  should  be  aware  that  only  8  panels  are  not 
sufficient  to  describe  the  geometry  in  most  of  the  cases. 
Basically  the  achieved  accuracy  depends  on  both  the  shape  of 
the  body  and  the  panel  configuration  (number  of  panels  and 
local  widths).  A  closer  spacing  is  advisable  in  regions  where 
severe  changes  of  the  pressure  distribution  are  expected 
(e.g.,  leading  edge). 


Flow  past  a  pair  of  circular  cylinders 


Two  circular  cylinders  are  arranged  side-by-side  in  a  uni¬ 
form  stream.  The  surface  of  each  cylinder  is  replaced  by  50 
panels  of  equal  width  (see  Figure  3.7(a)).  The  computed 
velocity  distribution  on  one  of  the  cylinders  is  shown  in 
Figure  3.8(b).  The  reader  shall  pay  some  attention  to  a 
comparison  between  the  flow  past  one  and  the  flow  past  a  pair 
of  cylinders.  Obviously  the  maximum  velocity  is  increased 
by  the  existence  of  a  second  cylinder.  The  closer  the  two 
cylinders  are  arranged,  the  higher  the  maximum  velocity. 

While  the  stagnation  points  in  a  single  cylinder  flow  are 
located  at  the  farthest  down  and  upstream  points  of  the  cylinder, 
the  disturbance  of  a  second  cylinder  causes  the  stagnation 
points  to  move  towards  the  other  cylinder.  The  streamline 
picture,  given  in  Figure  3.8,  should  provide  a  deeper  under¬ 
standing  of  this  kind  of  flow. 

Flow  past  two  element  airfoil 

The  main  goal  of  leading  and  trailing  edge  devices  is  to 
obtain  a  higher  lift  coefficient.  We  will  investigate  the 
effect  of  a  single  slotted  flap  on  the  pressure  distribution 
of  the  main  airfoil. 

The  pressure  distributions  of  a  single  airfoil  and  of  an 
airfoil-flap  combination  are  compared  in  Figure  3.9.  The 
coordinates  of  both  main  airfoil  (a  NACA  4412)  and  flap  are 
listed  in  Section  F  (sample  input  data)  .  The  results  indicate 
that  lift  increases  more  than  50%  by  using  a  slotted,  21.5 
degree  deflected  flap. 


I 


Figure  3.10).  Since  an  axis  of  symmetry  must  be  impermeable 
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on  the  upper  surface  and  increased  on  the  lower  surface.  In 
the  particular  case  the  overall  lift  gain  is  about  15%  of  the 
lift  in  free  air,  but  we  might  not  always  expect  a  lift  gain. 
The  actual  balance  between  lift  reduction  on  the  upper  and 
lift  increase  on  the  lower  surface  depends  on  both  distance 
from  the  ground  and  angle  of  incidence.  Figure  3.12  confirms 
that  there  are  cases  with  less  lift  than  in  free  air.  High 
angles  of  attack  and  moderate  distances  from  ground  are  sus¬ 
ceptible  constellations  to  lift  loss. 

F.  I/O — DESCRIPTION  AND  LISTING  OF  THE  PROGRAM  "PANEL" 

This  program  calculates  non-lifting  and  lifting  potential 
flow  past  one  or  two  bodies.  Any  2-dimensional  shape  and 
any  angle  of  attack,  which  do  not  cause  flow  separation, 
are  acceptable. 

Input  data 

The  data  must  be  arranged  in  the  following  order: 

(1)  Header  card; 

(2)  Coordinates  of  first  body  cards  (variable  number  of 
cards) ; 

(3)  Second  body  control  card; 

(4)  Coordinates  of  second  body  cards  (variable  number  of 
cards) . 

Items  3  and  4  are  used  only  for  the  2-body  case.  The  actual 
instructions  are  as  follows. 

Header  card 


1-10:  Number  of  bodies  (integer) 

11-20:  Number  of  points  of  the  first  body  (integer) 

21-30:  Angle  of  attack  in  degrees  (real) 


PANEL  is  the  number  of  the  panel; 

X  and  Y  are  the  coordinates  of  control  points  (not 
boundary  points); 

V  denotes  the  relative  velocity  (V/V^) ;  and 
C  denotes  the  pressure  coefficient. 
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NQNLIFTING  SOLUTION 


PANEL 

1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 

19 

20 
21 
22 

23 

24 

25 

26 

27 

28 

29 

30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 

41 

42 


XM 

YM 

0.97500 

-0.00080 

0.92500 

-0.00190 

0.85000 

-0.00305 

0.75000 

-0.00520 

0.65000 

-0.00825 

0.55000 

-0.01200 

0.45000 

-0.01600 

0.35000 

-0.02030 

0.27500 

-0.02380 

0.22500 

-0.02620 

0.17500 

-0.02810 

0.12500 

-0.02870 

0.08750 

-0.02800 

0.06250 

-0.02615 

0.03750 

-0.02220 

0.01875 

-0.01690 

0.00625 

-0.00715 

0.00625 

0.01200 

0.01875 

0.02895 

0.03750 

0.04060 

0.06250 

0.05245 

0.08750 

0.06175 

0.12500 

0.07240 

0.17500 

0.08345 

0.22500 

0.09105 

0.27500 

0.09585 

0.35000 

0.09780 

0.45000 

0.09495 

0.55000 

0.08665 

0.65000 

0.07415 

0.75000 

0.05790 

0.85000 

0.03800 

0.92500 

0.02090 

0.97500 

0.00735 

1.22500 

-0.14500 

1.17500 

-0.13250 

1.10000 

-0.10250 

1.02500 

-0.06500 

1.02500 

-0.05500 

1.10000 

-0.07500 

1.17500 

-0.10250 

1.22500 

-0.13250 

V 

CP 

-3.05132 

-8.31055 

-1.71788 

-1.95110 

-1.43382 

-1.05583 

-1.20962 

-0.46318 

-1.11906 

-0.25229 

-1.07063 

-0.14625 

-1.03543 

-0.07210 

-1.01792 

-0.03616 

-0.97360 

0.05210 

-0.98044 

0.03874 

-0.97498 

0.04942 

-0.95981 

0.07876 

-0.92944 

0.13613 

-0.90437 

0.18212 

-0.81613 

0.33393 

-0.66716 

0.55489 

-0.06604 

0.99564 

1.20372 

-0.44893 

1.51674 

-1.30050 

1.46366 

-1.14230 

1.45104 

-1.10550 

1.44313 

-1.08263 

1.40582 

-0.97632 

1.38729 

-0.92458 

1.36530 

-0.86405 

1.34794 

-0.81694 

1.24546 

-0.55117 

1.15851 

-0.34214 

1.07036 

-0.14566 

0.96331 

0.07203 

0.82090 

0.32612 

0.55160 

0.69573 

0.37152 

0.86197 

-0.86208 

0.25681 

-2.69341 

-6.25447 

-1.67415 

-1.80278 

-0.59324 

0  .64806 

1.34805 

-0.81725 

3.04627 

-8.27976 

1.49617 

-1.23852 

0.51932 

0.73031 

-1.15598 

-0.33629 

1. 
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LIFTING  SOLUTION 


PANEL 

XM 

YM 

V 

CP 

1 

0.97500 

-0.00080 

-1.09525 

-0.19958 

2 

0.92500 

-0.00190 

-0.80975 

0.34430 

3 

0.85000 

-0.00305 

-0.71356 

0.49083 

4 

0.75000 

-0.00520 

-0.65850 

0.56638 

5 

0.65000 

-0.00825 

-0.63122 

0.60157 

6 

0.55000 

-0.01200 

-0.60472 

0.63432 

7 

0.45000 

-0.01600 

-0.56660 

0.67897 

8 

0.35000 

-0.02030 

-0.52361 

0.72583 

9 

0.27500 

-0.02380 

-0.47182 

0.77738 

10 

0.22500 

-0.02620 

-0.43188 

0.81348 

11 

0.17500 

-0.02810 

-0.35514 

0.87387 

12 

0.12500 

-0.02870 

-0.22992 

0.94714 

13 

0.08750 

-0.02800 

-0.07829 

0.99387 

14 

0.06250 

-0.02615 

0.09794 

0.99041 

15 

0.03750 

-0.02220 

0.47846 

0.77108 

16 

0.01875 

-0.01690 

1.09475 

-0.19848 

17 

0.00625 

-0.00715 

2.42874 

-4.89880 

18 

0.00625 

0.01200 

3.46087 

-10.97764 

19 

0.01875 

0.02895 

3.13743 

-8.84344 

20 

0.03750 

0.04060 

2.71609 

-6.37716 

21 

0.06250 

0.05245 

2.48044 

-5.15257 

22 

0.08750 

0.06175 

2.34667 

-4.50686 

23 

0.12500 

0.07240 

2.20878 

-3.87871 

24 

0.17500 

0.08345 

2.08988 

-3.36759 

25 

0.22500 

0.09105 

2.00707 

-3.02833 

26 

0.27500 

0.09585 

1.94596 

-2.78674 

27 

0.35000 

0.09780 

1.81808 

-2.30540 

28 

0.45000 

0.09495 

1.69455 

-1.87149 

29 

0.55000 

0.08665 

1.60093 

-1.56299 

30 

0.65000 

0.07415 

1.51477 

-1.29451 

31 

0.75000 

0.05790 

1.43499 

-1.05920 

32 

0.85000 

0.03800 

1.33681 

-0.78706 

33 

0.92500 

0.02090 

1.28943 

-0.66263 

34 

0.97500 

0.00735 

1.09525 

-0.19956 

35 

1.22500 

-0.14500 

-0.86324 

0.25482 

36 

1.17500 

-0.13250 

-0.79007 

0.37579 

37 

1.10000 

-0.10250 

-0.43515 

0.81064 

38 

1.02500 

-0.06500 

0.41219 

0.83010 

39 

1.02500 

-0.05500 

1.91951 

-2.68453 

40 

1.10000 

-0.07500 

1.49954 

-1.24862 

41 

1.17500 

-0.10250 

1.23353 

-0.52159 

42 

1.22500 

-0.13250 

0.86324 

0.25482 
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Program  listing 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  THIS  PROGRAM  CALCULATES  2-D  POTENTIAL  FLOW  PAST  1  OR  2-BODIES  C 

C  AT  ANY  ANGLE  OF  ATTACK.  C 

C  C 

C  WRITTEN  BY  :  CAPT  LEE, HEE  WOO  C 

C  DATE  :  NOV. 28  1985  C 

C  C 

C  NOTE  i  MAXIMUM  NUMBER  OF  PANELS  =  200 
C 


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DIMENSION  Z(200,200),XBC200),YB(200),BEC200),TH(200),V(200) 
DIMENSION  WKAREA( 65000 ), XMC 200 ),YM( 200 ),VVVC 200 ),VCC 200) 

DIMENSION  CPC  200 ), VVC 200 ),YC 200, 200), AC  200), VI (200), SI (ZOO) 
DIMENSION  VC1(200),VC2(200),V2(200),VT(200) 

C 

C  -  READ  INPUT  DATA  (FOR  FIRST  BODY)  — 

C 

READC  4,1)  NB,NN, AN 

I  FORMAT (2110, F10. 5) 

AN  =  ANX3. 141592/180. 

DO  10  I  =  l.NN 

READ(4,11)XB(I),YB(I) 

II  FORMAT  C  2F10 . 5 ) 

10  CONTINUE 

C 

C  -  CALCULATE  MID-POINTS  OF  PANELS  AND  ANGLES  THETA  - 

C 

N  =  NN-1 
DO  12  1=1, N 
K  =  1+1 

XMC I )  =  (XBC I )+XB(K) )/2. 

YM( I )  =  CYB(I)+YB(K))/2. 

THH  =  CYBCK)-YBCI))/(XB(K)-XB(I)) 

THC I ) =  ATANCTHH) 

IFCXB(K) .LT.XB(I))  TH( I ) =TH( I )+3 . 141592 
12  CONTINUE 
C 

C  -  CALCULATE  PANEL  LENGTHS  - 

C 

DO  13  I  =  1 , N 
K  =  1  +  1 

A( I )  =  SQRT((XB(K)-XB(I))*X2+(YB(K)-YB(I))*x2) 

13  CONTINUE 

IF(NB.NE.l)  GO  TO  600 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  NON-LIFTING  PART  (1-BODY)  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
DO  14  I  =  1 ,  N 
C 

c  —  CALCULATE  ANGLE  BETA  AND  NORMAL  COMPONENT  OF  FREE  STREAM 
C  VELOCITY  — 

C 

BEC I ) =  TH(I)+3. 141592/2. 

VCD  =  -(COSCBE(I) )XC0SCAN)+SIN(BE(I))XSINCAN)) 

DO  15  J  =  1,N 


UFrr?ODUCED  AT  GOV ERHMfil  I 


I'v'.'v  i,ti  pr^.’w?yv*n.''-_i'v'  k-v 


C  —  CALCULATE  INFLUENCE  COEFFICIENTS  OF  NORMAL  VELOCITY  - 

C 

IF(I.EQ.J)  GO  TO  15 

B  =  (XM(I)-XB(J))XC0S(BECI))+(YM(I)-YB(J))XSINCBECI)) 

C  =  COS(TH(J))xCOS(BECI))+SIN(TH( J ) )XSINC BE( I ) ) 

E  =  2 . XCOSCTHC  J ))X(XM(I)-XB(J))+2.XSIN(TH(J))X(YM( I )-YB( J ) ) 
F  =  CXM(I)-XB(J))XX2+(YM(I)-YB(J))XX2 
H  =  SORT ( 4 . XF-EXX2) 

Z(I,J)  =  TEG(B,C,E,F,H,A(J)) 

15  CONTINUE 
ZCI,I)=3. 141592 

14  CONTINUE 
C 

C  -  SOLVE  SET  OF  LINEAR  EQUATIONS  FOR  SOURCE  STRENGTHS  - 

C 

IDGT  =  0 

CALL  LEQT2F  CZ, 1,N, 200, V, IDGT, WKAREA,  IER) 

C 

C  -  CALCULATE  INFLUENCE  COEFFICIENTS  OF  TANGENTIAL  VELOCITY  - 

C 

DO  16  I  =  1 ,  N 

DO  17  J  =  1 .  N 

IFCI.EQ.J)  GO  TO  17 

B  =  (XM(I)-XBCJ))XCOSCTHCI))+CYMCI)-YBCJ))XSIN(THCI)) 

C  =  COS(TH(J))xCOS(TH(I))+SIN(TH(J))XSIN(THCI)) 

E  =  2 . XCOSC  TH( J ))X(XM(I)-XB(J))+2.XSIN(TH(J))X(YM( I )-YB( J ) ) 
F  =  (XM(I)-XBCJ))x*2+CYM(I)-YBCJ))xx2 
H  =  SQRT ( 4 . XF-EXX2) 

Y( I , J )  =  TEGCB,C,E,F,H,ACJ)) 

17  CONTINUE 

Y(I,I)=0.0 

16  CONTINUE 
C 

C  -  CALCULATE  TOTAL  VELOCITY  AND  CP  AT  MIDPOINTS  OF  EACH  PANEL  - 

C 

WRITEC8.95) 

95  F0RMATC///.25X, ’NONLIFTING  SOLUTION*,// 

xiOX, ’PANEL', 5X, *XM*,8X, *YM*,11X, *V’,10X, *CP’) 

DO  18  I  =  1 , N 
S  =  0. 

DO  19  J  =  1 , N 

S  =  S+V(J)XY(I, J) 

19  CONTINUE 

VVCI)  =  COS(TH(I))XCOS(AN)+SIN(THCI))XSIN(AN)+S 
CPCI)  =  l.-VV(I)xx2 
WRITEC  8 ,  93 )  I,XM(I),YM(I)»VV(I),CP(I) 

93  FORMAT (10X,I3,3X,2F10.5,3X,2F10.5) 

18  CONTINUE 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c  c 

C  LIFTING  PART  (1-BODY)  C 

C  C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  —  CALCULATE  SOURCE  STRENGTHS  DUE  TO  CIRCULATORY  FLOW  OF  UNIT 


STRENGTH  - 


DO  20  I  =  1,N 

$  =  0. 

DO  21  J  =  1 ,  N 

S  =  S+Y ( I , J  ) 

21  CONTINUE 

V1(I)=S 

20  CONTINUE 
IDGT  =  0 

CALL  LEQT2F  CZ, 1 , N, 200 , VI , IDGT, WKAREA, IER) 


'  •  «  •  »  •  •  •»•»*•*•*»  4**  ' 


-  CALCULATE  DISTURBANCE  VELOCITY  DUE  TO  CIRCULATORY  FLOW  OF  UNIT 

STRENGTH  - 

DO  22  I  =  1,N 
S  =  0. 

DO  23  J  =  1,N 

S  =  S+YCI, J)*V1(J)+ZCI,  J) 

23  CONTINUE 

VC(I)=S 

22  CONTINUE 

-  CALCULATE  VORTEX  STRENGTHS  BY  KUTTA  CONDITION  - 

SI  =  -( VV( 1 )+VV(N) )/( VC(1 )+VC(N) ) 

-  CALCULATE  TOTAL  VELOCITY  AND  CP  AT  MIDPOINTS  OF  EACH  PANEL  - 

WRITE(8,96) 

96  F0RMAT(///,27X, 'LIFTING  SOLUTION',// 

K10X, 'PANEL ',5X, 'XM',8X, *YM',11X, 'V* , 10X, 'CP» ) 

DO  24  I  =1 ,  N 

vvvci)  =  vvcn+si*vca) 

CP(I)  =  1 . -VVVCI )**2 

WRITE(8,93)I,  XM(I),YM(I), VVVCI ) ,CP( I ) 

24  CONTINUE 

25  FORMAT ( 3F10 .5) 

GO  TO  700 

-  READ  INPUT  DATA  (FOR  SECOND  BODY)  - 

600  READ( 4,31)  MM 

31  FORMAT ( 110  ) 

DO  32  I  =  NN+1 , NN+MM 

READ(4,11)XB(I),YBCI) 

32  CONTINUE 

M  =  NN+MM-2 

-  CALCULATE  MID-POINTS  OF  PANELS  AND  ANGLES  THETA  — 

DO  33  I=N+1,M 
K  =  1  +  1 

XMC I )  =  (XB(K)+XB(K+l))/2. 

YM( I )  =  (YB(K)+YB(K+l))/2. 

THH  =  (YB(K+1)-YBCK))/(XB(K+1)-XB(K)) 

TH( I ) =  ATAN(THH) 

IF(XB(K+1) .LT.XB(K))  TH( I ) =TH( I )+3 . 141592 

33  CONTINUE 

-  CALCULATE  PANEL  LENGTHS  - 

DO  34  I  =  N+l , M 
K  =  1  +  1 

A( I )  =  SQRT((XB(K+l)-XB(K>)x*2+CYBCK+l)-YBCK))**2) 

34  CONTINUE 

DO  35  I  =  N+l , M+l 

XBCI)  =  XB( 1+1 ) 

YBC I )  =  YBCI+l) 

35  CONTINUE 


ooo  ooo  ooo  ooo  noon  oooo 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

NON-LIFTING  PART  (2-BODIES) 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

DO  36  I  =  1 » M 


-  CALCULATE  ANGLE  BETA  AND  NORMAL  COMPONENT  OF  FREE  STREAM 

VELOCITY  — 


BE(I)=  TH(I)+3. 141592/2. 

V(I)  =  -CCOS(BECI))xCOS(AN)+SINCBECI))xSIN(AN)) 
-  CALCULATE  INFLUENCE  COEFFICIENTS  OF  NORMAL  VELOCITY  - 


DO  37  J  =  1,M 

IFCI.EQ.J)  GO  TO  37 

B  =  CXMCI)-XBCJ))XCOSCBECI))+CYMCI)-YB(J))XSIN(BECI)) 

C  =  C0SCTHCJ))XC0SCBECI))+SIN(THCJ))XSINCBECI)) 

E  =  2.*C0S(TH(J))X(XM(I)-XB(J))+2.XSIN(TH(J) ) X( YMC I)-YB(J)) 
F  =  (XM(I)-XB(J) )XX2+(YM(I)-YB(J))**2 
H  =  SQRT(4.xF-Ex*2) 

ZCI,J)  =  TEG(B,C,E,F,H,ACJ)) 

7  CONTINUE 

Z(I,I)=3. 141592 
6  CONTINUE 

—  SOLVE  SET  OF  LINEAR  EQUATIONS  FOR  SOURCE  STRENGTHS  - 


IDGT  =  0. 

CALL  LEQT2F  (Z, 1 , M, 200 , V, IDGT, WKAREA , IER) 

-  CALCULATE  INFLUENCE- COEFFICIENTS  OF  TANGENTIAL  VELOCITY  - 

DO  38  I  =  LH 

DO  39  J  =  1,M 

IFCI.EQ.J)  GO  TO  39 

B  =  <  XMC I )-XB( J ) )XCOS(TH( I ) )  +  (YM( I )-YB( J ) )XSIN(TH(I)) 

C  =  C0S(TH(J))XC0S(TH(I))+SIN(TH(J))XSIN(THCI)) 

E  =  2 . XC05 (TH( J ) )X(XM( I )-XB( J ) )+2 . XSIN( TH( J ) ) x( YMC I)-YB(J)) 
F  =  CXMCI)-XB(J)) XX2+C YMC I)-YBCJ))x*2 
H  =  SQRTC4.XF-EXX2) 

Y( I , J  )  =  TEG(B,C,E,F,H,ACJ)) 

39  CONTINUE 

YCI,I)=0.0 

33  CONTINUE 
WRITEC8, 95) 

-  CALCULATE  TOTAL  VELOCITY  AND  CP  AT  MIDPOINTS  OF  EACH  PANEL 

DO  40  I  =  1,M 
S  =  0. 

DO  41  J  =  1 ,  M 

S  =  S+VC  J )XY( I , J  ) 

41  CONTINUE 

VVCI)  =  COSCTHC I ) )XCOSCAN)+SINCTHC I ) )XSINC AN)+S 
CP(I)  =  1.-VVCDXX2 
WRITEC8 , 93 )  I , XMC I) , YMC I ) , VVC I ) , CPC  I ) 

40  CONTINUE 
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ooo 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 
c 

C  LIFTING  PART  (2-BODIES) 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C  —  CALCULATE  SOURCE  STRENGTHS  DUE  TO  CIRCULATORY  FLOW  OF  UNIT 

C  STRENGTH  PAST  FIRST  BODY  - 

C 

DO  42  I  =  1,M 
S  =  0. 

DO  43  J  =  1,M 
G  =  1. 

IF(J.GT.N)  G=0 . 

S  =  S+Y( I , J ) *G 

43  CONTINUE 

V1(I)=S 

42  CONTINUE 
IDGT  =0. 

CALL  LEQT2F  (Z, 1 , M, 200 , VI , I DGT , WKAREA ,  IER) 

C 

C  —  CALCULATE  DISTURBANCE  VELOCITY  DUE  TO  CIRCULATORY  FLOW  OF  UNIT 

C  STRENGTH  PAST  FIRST  BODY  - 

C 

DO  44  I  =  1,M 
S  =  0. 

DO  45  J  =  I'M 
G  =  1. 

IFCJ.GT.N)  G=0 . 

S  =  S+YCI, J)XV1CJ)+Z(I,J)XG 

45  CONTINUE 

VC1CI)=S 

44  CONTINUE 
C 

C  ---  CALCULATE  SOURCE  STRENGTHS  DUE  TO  CIRCULATORY  FLOW  OF  UNIT 
C  STRENGTH  PAST  SECOND  BODY  — 


DO  46  I  =  1,M 
S  =  0. 

DO  47  J  =  1,M 
G  =  1. 

IFCJ.LE.N)  G=0. 

S  =  S+Y(I,J)XG 

47  CONTINUE 

V2C I ) =S 

46  CONTINUE 
I DGT=0 . 

CALL  LEQT2F  (Z,  1 , M,  200 , V2, IDGT, WKAREA , I ER) 

C 

C  —  CALCULATE  DISTURBANCE  VELOCITY  DUE  TO  CIRCULATORY  FLOW  OF  UNIT 
C  STRENGTH  PAST  SECOND  BODY  - 


DO  48  I  =  l.M 
S  =  0. 

DO  49  J  =  l.M 
G  =  1. 

IFCJ.LE.N)  G=0 . 

S  =  S+Y(I, J)XV2(J)+ZCI, J)XG 

49  CONTINUE 

VC2( I  )=S 

48  CONTINUE 


« 
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ooo 


ooo  nno  non 


-  CALCULATE  VORTEX  STRENGTHS  BY  KUTTA  CONDITION  — 

XI  =  VC1(1)+VC1(N) 

X2  =  VC2C 1 )+VC2( N ) 

X3  =  -Wd)-VV(N) 

X4  =  VC1(N+1)+VC1(M) 

X5  =  VC2(N+1 )+VC2(M) 

X6  =  -VV(N+1)-VVCM) 

SI2  =  CX3xX4-XlxX6)/(X4xX2-Xl*X5) 

SI1  =  CX3-X2XSI2)/X1 

-  CALCULATE  TOTAL  VELOCITY  AND  CP  AT  MIDPOINTS  OF  EACH  PANEL  — 

WRITEC8  > 96 ) 

DO  50  I  =1,M 

VT(  I )  =  VCl(nxSIl+VC2CI)XSI2 
vvvci)  =  wcn+VT(i) 

CP(I)  =  1.-VVV(I)XX2 
WRITE(8,93)I,XM(n,YM(I),VVV(I),CPCI) 

50  CONTINUE 
00  WRITEC6  > 97  ) 

97  FORMATC1X, 'COMPUTATION  COMPLETED') 

STOP 

END 

-  THIS  FUNTION  EVALUATES  THE  INTEGRALS  (INFLUENCE  COEFFICIENTS) 

FUNCTION  TEG(B,C,P,Q,R,S) 

TERM1  =  AL0G((SXX2-PXS+Q)/Q) 

TERM2  =  ATAN((2.XS-P)/R)-ATAN(-P/R) 

TEG  =  -CxTERMl/2.+(2.XB-CxP)XTERM2/R 

RETURN 

END 
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IV.  BOX  METHOD 


A.  INTRODUCTION 

The  thin  shear  layer  equations  are  more  complicated  than 
Laplace's  equation  because  they  are  nonlinear.  This  chapter 
presents  the  box-method,  which  can  be  applied  to  the  solution 
of  the  thin  shear  layer  equations.  The  box  method  was  intro¬ 
duced  by  Keller  in  1970  [Ref.  4]. 

One  of  the  basic  ideas  of  the  box  method  is  to  write  the 
governing  system  of  equations  in  the  form  of  a  first-order 
system.  This  system  is  solved  by  finite-difference  approxima¬ 
tions  and  Newton's  method  is  applied  to  solve  the  equations. 
Finally,  the  resulting  linear  system  is  solved  by  the  block- 
elimination  method. 

B.  FALKNER-SKAN  TRANSFORMATION 

The  thin  shear  layer  equations  for  incompressible  laminar 
flow  take  the  form 


3u  3v 
3x  3y 


0 


(4.1) 


u 


v 


3u 

3y 


1  dP 
p  dx 


(4.2) 


Boundary  conditions  are  prescribed  at  the  surface 


and  at  the  edge  of  the  boundary  layer 


y  -*■  *  u  =  U  (x) 

A  e 


(4.3) 


It  is  convenient  to  reformulate  the  equations  using  the 
streamfunction  and  the  similarity  concept.  Therefore  the 
Falkner-Skan  transformation  is  introduced. 


n 


U  1/2 

,  e> 

W  y 


Re 


1/2 


(4.4) 


ip  (x,y) 


(U 


vx)  1//2f  (x,n) 


(4.5) 


Substituting  Eqs.  (4.4)  and  (4.5)  into  Eq.  (4.2),  we  get  the 
transformed  momentum  equation  for  2-D  laminar  flows. 


f." 


+  52Jiff"  +m  [1-  (f  '  )  2  ] 


=  x  ( f 


9x 


9x 


(4.6) 


where 


m 


U  dx 
e 

(dimensionless  pressure-gradient) 


with  the  boundary  conditions 


UFPfJOUUCED  AT  GO  V  F  P I  J»1  r r 1  l 


If  f  is  a  function  of  n  only,  the  right-hand  terms  of  Eq. 

(4.6)  will  be  zero.  Then  this  will  be  a  third-order  ordinary 
differential  equation  whose  solution  is  well-known  as  a 
"similar  flow."  But,  if  f  is  a  function  of  n  and  x  (non¬ 
similar  flows),  we  need  a  numerical  procedure,  such  as  the 
box  method. 

C.  NUMERICAL  FORMULATION  (BOX  METHOD) 

First  of  all,  the  coordinates  (x,y)  of  a  given  geometry 
must  be  transformed  to  coordinates  (£,n)  to  apply  the  Box 
method  (see  Figure  4.1). 


Stagnation  Point 

T  • 


Ai  rfoi 1 


_ 

^77rr 

(/  BOX'-- 

W//.  ■■  ■ 

j 

i 

&  **  . 


£  § 


Figure  4.1.  Transformed  Coordinates  of  Upper  Surface 
Airfoil  and  Net  Rectangle  for  Difference 
Approximations 


-  •  -A-  .  •  o  .  •  .  •  .  ■  - 


*  *> 


The  boundary  layer  thickness  in  transformed  coordinates  is 


nearly  independent  of  the  streamwise  distance  and  can  be 
represented  by  a  fixed  number  of  profile  points  at  fixed 
spacing. 

One  of  the  basic  ideas  of  the  Box  method  is  to  write  the 
governing  system  of  equations  in  the  form  of  a  first  order 
system.  We  write  Eq.  (4.6)  in  terms  of  a  first-order  system 
of  PDE's 


f  '  =  u  ( £  ,  n  ) 

u'  =  v(£,n) 

(bv)  '  +  (2ii)fv  +  m  (1-u2)  =  £(u  -  v  ||) 

where  a  prime  denotes  differentiation  with  respect  to  n 
and  b  =  1  +  v^/v 

with  the  boundary  condition 

f(£,0)  =  0,  u  ( £ ,  0 )  =  0,  u(£,noo)  =  1 

We  denote  the  net  points  shown  in  Figure  4.1  as 


II 

o 

0 

F  . 

l 

=  £ .  i  +  k . 

l-l  l 

i  =  1,2,  . 

.  .  » I 

S3 

O 

II 

0 

nD 

=  n  •  ••  +  h . 

3-1 

j  =  1,2,. 

•  »  /  (J 

(4.8a) 

(4.8b) 

(4.8c) 


(4.9) 
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And  we  can  introduce  the  following'  approximations : 


A)  Coordinates  of  midpoints  (£  ^,o  ^)  and  net  functions 


(g  stands  for  f,  u  or  v) 


1-o  3“o 


5  1  =  n.  1 


1-2 


I(nj  +nj-l) 


(4.10a) 


.  1 

X~J  -  1/  i  ^  i-l>  i  -  1/  i  ,  i  \ 

g-i  =  j(gj  +gj  }  9-  I  =  y  9h  g^-i} 

j"2 


2^j 


(4.10b) 


where  [  ] ^  means  the  quantities  (f  or  u  or  v)  at  point  (5^,n-) 


B)  Finite-difference  approximation 

From  Eqs .  (4.8a)  and  (4.8b),  the  centered-dif f erence 

derivatives  are 


f1  -f1  , 

.2  ...jr1. 


-  u 


.  1 
3  "2 


(4.11a) 


l  i 
u  .  -  u  .  , 

_J _ 2Z1 


=  v 


1 
3~2 


(4.11b) 


After  introducing  these  approximations  into  Eq.  (4.8)  and 

rearranging  (the  known  quantities  are  moved  to  the  right  hand 

side),  we  get  the  equation  (4.12c)  which  is  centered  about  the 

point  (£  lfn  x) .  This  represents  the  relationship  of  quan- 
i-2  ^~2 

tities  between  the  points  of  the  box. 

h .  . 

f1  -  f1  .  - ^(u1  +  u  ,)  =  0  (4.12a) 


0 


(4.12b) 


uj  -  Vi  -  -*<v; +  vi>  - 


i  i  Li 


.  2 .  i 


where 


a  = 


,i-l 


i-1 


+  a  (v 


' .  1 
1~2 


-L  3 

i  .  i 

3-2 

z  .X 

3-2 

i-l_i 

-i-1  i  . 
f.  iv.  i> 

_  i — 1 

.  lf.  1  - 

=  R.  1 

3  2  3  2 

j—  j-j 

m .  +1 

3-2 

p 
t— * 

ii 

—^2 —  +  a' 

a2  =  mi 

i-1 


\  1 
3“2 


{-L  x  +a[  (fv)  j-fu*)  x]  }  -in 


3  "7 


3-2  3-2 


(4.12c) 


J.  1 
3-2 


b.v.-b.  v ,  ,  ,  i  ^  i-1 

3-1  +  g^L(fv)  1+m[l-(u2)  ]} 

3  3-2 


3  ~2 


The  last  of  the  above  equations  is  non-linear.  Therefore  we 
introduce  Newton's  method  to  solve  this  system.  We  set 


gin+1) 


+  6g  . 
J 


(n) 


n  =  0,1,2,... 


(4.13) 


where  the  superscript  in  parentheses  is  the  iteration  counter 
with  initial  condition 


B 


Introducing  Eq.  (4.13)  to  Eq.  (4.12)  and  dropping  the  quadratic 
terms  in  g^,  we  get  (superscripts  i  and  n  are  dropped  for 
simplicity) 


sfi  ‘  6Ej-i  -  4(<uj  +  6uj-i»  *  <ri»j 


(4.15a) 


■j  •  Suj-i  *  4<6vj  +  - 


(4.15b) 


(VjSVj  +  (S2).6v..1+  (S3)j6f.  +  (S4)j6f._1 


+  (S5>jauj  +  (S6)j6uj-l  =  <r2> j 


where  all  terms  are  explained  in  Reference  5,  with  the 
boundary  conditions 


(4.15c) 


<5f  q  =  0  6uq  =  0  6Uj  =  0 


(4.16) 


D.  BLOCK  ELIMINATION  METHOD 

This  is  a  very  effective  way  to  solve  linear  difference 
equations,  discussed  by  Keller  in  1974  [Ref.  6].  We  write 


Eq.  (4.15)  in  a  matrix-vector  form 


A6  =  r 


where 


A  = 


^0  CC 


B. 


*1 


BJ-1  AJ-1  CJ-1 
B,  A, 


where 


B. 

3 


-1  -h ./2 
3 

(S4)j  (Vj 


(s2)j 


1  0 

0 
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0  1 

0 
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3 
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(s5)j 

<S1>J 
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3 

0 

-1 

-hj+l 

1 

(S 

0 

0 

0 
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1  <  j 


where 


I  is  the  identity  matrix 


I  = 


10  0 
0  10 
0  0  1 


From  Eq.  (4.18),  we  find  that 


Qo  -  Ao 


(4.19a) 


Qj 


B  . 
1 


j  1,2,.. .,J 


Aj~PjCj_^  j  —  1,2, ...,J 


(4.19b) 


(4.19c) 


I 

8 


Keller  showed  that  the  matrix  Pj  has  the  same  structure  as 
the  matrix  B ^ .  From  Eq.  (4.19),  we  derive  the  elements  of 
Pj  and  Qj . 
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o  <  j  <  J-l 
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V.  INTERACTION  METHOD 


A.  INTRODUCTION 

Interactive  methods  provide  a  special  coupling  between 
viscous  and  inviscid  flows.  They  are  intended  to  compute 
flows  including  separation.  Thus  these  methods  may  be  regarded 
as  an  alternative  to  the  Navier-Stokes  solvers,  which  are 
hardly  engineering  tools  because  of  their  huge  computer  time 
and  storage  requirements. 

The  classical  method  to  compute  viscous  flows  past  airfoils 
proceeds  as  follows:  The  velocity  distribution  is  computed 
by  any  appropriate  inviscid  flow  solver.  The  output  of  the 
inviscid  flow  solver  becomes  the  input  of  the  viscous  flow 
solver.  Solving  for  viscous  flow  consists  of  integrating  the 
boundary  layer  equations.  Provided  that  the  flow  remains 
attached  this  procedure  allows  a  reliable  prediction  of  lift 
and  drag.  Information  is  transferred  only  once  from  inviscid 
to  viscous  regions.  However,  many  flows  cannot  be  modelled  by 
one-time  information  transfers. 

Separation  bubbles  and  separated  flows  especially  require 
a  close  coupling  between  viscous  and  inviscid  regions.  Inter¬ 
action  schemes  provide  a  better  exchange  of  information  between 
these  two  regions. 

The  elements  of  interaction  schemes  are:  direct  or  inverse 
inviscid  flow  solver  and  direct  or  inverse  viscous  flow  solver. 
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The  direct  boundary  layer  method  has  the  disadvantage  that 
the  boundary  layer  equations  are  singular  at  the  point  of 
separation.  However,  if  the  external  velocity  is  computed 
by  prescribing  a  displacement  thickness  (known  as  the  inverse 
boundary  layer  method) ,  they  can  be  integrated  through  the 
point  of  separation. 

The  next  problem  associated  with  the  regions  of  reversed 
flow  is  numerical  instability,  because  downstream  marching 
procedures  cannot  be  applied  in  regions  of  reversed  flow. 

The  most  common  approximation  to  get  this  instability  under 
control,  the  so-called  FLARE  approximation,  neglects  the 
momentum  transport  term  u  9u/3x  in  regions  of  reversed  flow  as 
long  as  this  region  is  small.  However,  as  the  size  of  this 
region  increases,  the  FLARE  approximation  becomes  inaccurate. 

One  of  the  procedures  which  can  be  taken  into  account  is 
called  the  DUIT  (Downstream,  Upstream  Iteration) .  It  consists 
of  a  sequence  of  alternating  up  and  downstream  sweeps. 

There  are  several  types  of  recently  developed  interaction 
models.  All  procedures  have  to  solve  both  the  inviscid  (Laplace 


equation)  and  viscous  flow,  whose  equations  can  be  written 
according  to 
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where 


b  =  l  =  constant  in  laminar  flow 
b  =  1  +  vt/v  in  turbulent  flow 

Four  interaction  models  can  be  distinguished:  Direct,  Inverse 
Semi-inverse,  and  Simultaneous  interaction  methods  which  are 
subject  to  different  boundary  conditions. 

The  most  advanced  interaction  scheme  is  the  simultaneous 
interaction.  We  call  it  the  "strong  interaction"  (direct  and 
inverse  interactions  guarantee  weak  coupling  only) .  Examples 
in  Section  V.D  are  computed  by  the  Cebeci  program  using  this 
method.  Good  agreement  is  obtained  between  the  results  of 
interaction  schemes  and  experimental  results. 

B.  FOUNDATION  OF  THE  INTERACTION  SCHEMES 
1 .  Direct  Interaction  Scheme 

The  direct  interaction  model  is  composed  of  a  direct 


inviscid  and  a  direct  viscous  flow  solver  (see  Figure  5.1a). 
The  usual  seauence  is: 
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Figure  5.1.  Global  Organization  of  Interaction  Methods 
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(1)  Calculate  the  external  velocity  distribution  by 
inviscid  flow  computations. 

(2)  Calculate  the  displacement  thickness,  6*,  by  viscous 
flow  computations  using  the  external  velocity  as 
boundary  condition. 

(3)  Compute  an  updated  shape  of  the  displacement  body  and 
repeat  steps  1  and  2  until  the  results  converge. 

Unfortunately,  the  direct  boundary  layer  method  suffers 

numerical  breakdown  at  the  point  of  separation  (Goldstein 

singularity) .  Therefore  this  scheme  is  not  appropriate  to 

handle  airfoil  flows  with  separation. 

2 .  Inverse  Interaction  Scheme 


This  method  was  introduced  to  overcome  the  singularity 
problems  near  separation.  It  combines  an  inverse  inviscid  and 
an  inverse  viscous  flow  solver  (see  Figure  5.1b).  However, 
the  overall  procedure  suffers  from  very  slow  convergence. 

Due  to  this  shortcoming  one  shall  apply  this  inverse  scheme 
to  regions  of  separated  flow  only. 

3.  Semi-Inverse  Interaction  Scheme 


This  method  combines  a  direct  inviscid  flow  solver 
with  an  inverse  viscous  flow  solver  with  the  same  input  (dis¬ 
placement  thickness) .  This  leads  to  two  external  velocity 

distributions,  U  T(x)  and  U  ,.(x)  (see  Figure  5.1c).  Satis- 

el  ev  J 

factory  convergence  is  ensured  by  a  relaxation  formula,  which 
is  introduced  to  define  an  updated  displacement  thickness 
distribution . 
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where  u  is  a  relaxation  parameter.  The  numerical  weakness  of 
the  purely  inverse  scheme  is  improved  by  this  method,  but 
both  inviscid  and  viscous  regions  are  still  coupled  loosely. 

4 .  Simultaneous  Interaction  Scheme 

The  simultaneous  interaction  scheme  emphasizes  strong 
interaction  between  the  outer  inviscid  and  the  inner  viscous 
region.  The  external  velocity  Ug(x)  and  the  displacement 
thickness  S*(x)  are  treated  as  unknown  quantities.  An  addi¬ 
tional  relation  is  added,  the  so-called  interaction  law  which 
can  be  given  by  the  "blowing  velocity"  concept. 

The  equations  are  solved  by  the  inverse  method  with 
successive  sweeps  over  the  airfoil  surface  (see  Figure  5. Id). 
This  method  is  compatible  with  the  weak  interaction  scheme 
where  both  inviscid  and  viscous  regions  are  coupled  loosely. 

For  each  sweep,  the  external  velocity  for  the  boundary 
layer  equation  is  written  as 

U  (x)  -  U°(x)  +  6U  (x)  (5.4) 

e  e  e 


where 

U°(x)  is  the  inviscid  velocity; 

<5Ue(x)  is  the  perturbation  due  to  the  dis¬ 
placement  effect  of  a  boundary  layer. 

The  blowing  velocity  concept  is  introduced  to  get  the  pertur¬ 
bation  velocity  6U  by  the  interaction  law.  The  displacement 


effect  of  a  boundary  layer  can  be  modelled  by  ejecting  fluid 
at  the  airfoil's  surface  (see  Figure  5.2). 


Figure  5.2.  Blowing  Velocity  Concept 


A  properly  arranged  source  distribution  on  the  surface  dis¬ 
places  the  streamlines  away  from  the  surface  such  that  the 
virtual  displacement  body  becomes  a  streamline. 

Our  first. goal  is  to  determine  the  source  strengths  such 
that  the  tangential  flow  condition  on  the  displacement  body 
takes  the  form 


v (x , 6  * )  _  d3  * 

U  (x)  dx 

e 


To  achieve  this  goal  we  use  the  thin  airfoil  approximation: 

(1)  The  displacement  thickness  is  assumed  to  be  so 
small  that  u-velocity  components  do  not  vary  across 
the  layer. 

(2)  The  airfoil  in  this  connection  can  be  represented  by 
a  straight  line.  This  approximation  implies  that 
the  blowing  velocity  v(x,0)  equals  half  of  the 
source  strength. 
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where  (Ue6*)  is  defined  as  blowing  velocity.  Our  second 
goal  is  to  relate  the  perturbation  velocity,  6Ue  to  the  blowing 
velocity.  This  process  is  quite  similar  to  evaluating  tangen¬ 
tial  velocities  in  the  panel  method.  In  fact,  this  is  even 
simpler  because  of  the  straight  line  surface. 
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where  the  interaction  region  is  limited  to  a  finite  range 
xa  <_  x  _<  x^.  This  integral  is  referred  to  as  the  Hilbert 
integral.  Rewriting  Eq .  (5.4)/  we  finally  obtain  the  inter¬ 


action  law. 


The  numerical  implementation  of  Eq.  (5.8)  requires  a  discrete 
approximation  of  the  Hilbert  integral.  This  can  be  performed 
by  using  the  trapezoidal  rule. 

The  examples  in  Section  V.D  demonstrate  that  this  inter¬ 
action  method  can  give  reliable  results  for  flows  up  to  high 
angle  of  attack,  including  flows  with  bubbles  and  separation. 


C.  CONSIDERATION  OF  BOUNDARY  LAYER  TRANSITION  AND  OF 
TURBULENT  FLOW  MODELLING 

1 .  Transition 

One  of  the  most  important  parameters  to  predict  the 
drag  and  lift  of  an  airfoil  is  the  transition  point.  Boundary 
layer  transition  is  affected  by  many  parameters,  for  example, 
the  pressure  distribution  (major  parameter) ,  the  wall  roughness 
and  the  intensity  of  the  free  stream  turbulence,  etc.  Because 
of  this  fact,  the  theoretical  modeling  of  transition  is  very 
complicated  and  one  therefore  resorts  to  experimental 
information. 

In  the  Cebeci  program,  the  following  experimental 
correlation  formula  is  used,  which  was  given  by  Cebeci  and 

Smith  (1974)  as  a  relation  between  R-  and  Re  at  transition. 

0  x 
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where 


Re  =  U  x/v 
x  e 


and  0  is  the  momentum  thickness. 


2 .  Turbulent  Flow  Model 

Unlike  laminar  flows,  turbulent  flows  have  a  compli¬ 
cated  time-dependent  behavior.  It  is  too  difficult  to  deal 
with  the  instantaneous  properties.  Thus,  the  mean-flow 
properties  are  applied  in  turbulent  flow. 

The  most  common  mean-flow  models  are  the  "eddy- 
viscosity"  formula  which  are  based  on  thin  shear  laye’r 
assumptions . 


-  pu'v'  =  pvfc  7y  (5.10) 

where  is  related  empirically  to  the  mean  flow  velocity 
gradient  and  the  length  scale.  In  the  Cebeci  program,  vfc 
is  presented  by  the  algebraic  eddy-viscosity  formulation  of 
Ceb«!ci  and  Smith. 
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More  detailed  descriptions  are  listed  in  Reference  7. 

D.  EXAMPLES 

The  subsequent  examples  were  computed  using  a  program 
developed  by  Cebeci  and  coworkers  [Ref.  17] ,  on  the  NPS  IBM  370 
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1 .  Demonstration  of  the  Program  Capabilities 

The  velocity  profiles  on  both  upper  and  lower  surfaces, 
as  well  as  in  the  wake,  are  presented  in  Figure  5.3.  At 
this  angle  of  attack  (a  =  10°)  ,  transition  occurs  very  close 
to  the  leading  edge  on  the  upper  surface.  The  boundary  layer 
thickness  is  quite  thin  in  the  accelerated  flow  region  (right 
after  the  leading  edge) ,  but  it  grows  thicker  farther  down¬ 
stream  in  the  decelerated  flow  region  (near  the  trailing  edge) . 
Eventually,  we  find  a  small  reversed  flow  region  just  before 
the  trailing  edge  in  this  case.  The  wake  region  shows  the 
mixing  layer  which  decays  with  increasing  downstream  distance. 

Figure  5.4  demonstrates  how  lift,  drag  and  the  loca¬ 
tion  of  transition  depend  on  the  angle  of  attack.  The  skin 
friction  drag  is  dominant  at  low  angles  of  attack,  the  pressure 
drag  at  high  angles  of  attack  (sea  Figure  5.4b). 

Figure  5.5  shows  the  distributions  of  the  skin  friction 
coefficient,  displacement  thickness  and  shape  parameter  in 
dependence  of  Reynolds  number  and  angle  of  attack.  In  the 
attached  flow,  the  skin  friction  coefficient  decreases  along 
the  downstream  direction  until  the  point  of  transition 
(laminar  region) ,  but  increases  steeply  after  transition  and 
then  decreases  again  because  the  skin  friction  is  related  to 
the  slope  of  the  velocity  profile,  3u/9y,  at  the  surface.  At 
high  angles  of  attack,  transition  on  the  upper  surface  occurs 
close  to  the  leading  edge  and  the  negative  skin  friction 
coefficient  near  the  trailing  edge  indicates  separated  flow. 
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Also  the  displacement  thickness  and  the  shape parameter  increase 
drastically  in  regions  of  separated  flow  (see  Figure  5.5a). 

At  low  Reynolds  numbers,  laminar  flow  can  separate  at 
mid-chord  and  reattach  (we  call  this  phenomenon  a  "bubble"). 
Reattachment  can  often  occur  if  the  pressure  gradient  de¬ 
creases  rapidly  soon  after  separation,  so  that  a  strong  reversed 
flow  is  not  established.  Thus  the  shear  layer  reattaches  onto 
the  surface.  Accordingly,  the  displacement  thickness  and  the 
shape  parameter  increase  in  the  bubble  region  (see  Figure 
5.5b)  . 

2 .  Comparison  with  Experimental  Results 

The  comparison  of  pressure  distributions  is  shown  in 
Figure  5.6.  The  overall  lift  of  inviscid  calculations  devi¬ 
ates  20%  from  the  experimental  results.  The  interaction  method 
improves  the  accuracy,  but  the  computation  still  overestimates 
the  lift.  The  lift  coefficient  obtained  by  Cebeci ' s  program 
is  approximately  10%  larger  than  the  measured  one  (see  Reference 
8)  . 

Figure  5.7  demonstrates  that  the  accuracy  of  this  method 
drops  with  lower  Reynolds  numbers.  The  reasons  for  this  de¬ 
creased  accuracy  at  low  Reynolds  numbers  are  the  higher  proba¬ 
bility  to  separate  and  the  used  turbulence  model,  which  seems 
to  be  inappropriate  for  low  Reynolds  numbers.  Therefore,  low 
Reynolds  numbers  and  high  angles  of  attack  pose  severe  limi¬ 
tations  (see  Figure  5.7a  &  b) .  The  experimental  results  are 
taken  from  References  7  and  9. 
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. . :  Experimental 

:  Interaction  method 

Cp  and  Cd  Curves  for  (a)  FX  63-137,  Re  =  199487 
(b)  GA-1,  Re  =  1900000,  (c)  FX  63-137, Re  -  280000 
(  Source  of  experimental  results:  (a)Ref.9  (b)Ref.7  (c)Ref.  10 
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On  the  other  hand,  experimental  measurements  should 
not  be  expected  to  be  exact.  Turbulence  level  and  interference 
effects  influence  the  wind  tunnel  measurements.  Figure  5.7c 
shows  good  agreement  between  computed  and  experimental  results 
taken  from  Reference  10,  at  low  Reynolds  numbers. 

The  location  of  transition  and  separation  points  have 
an  important  influence  on  the  lift  and  drag  coefficients. 

Figure  5.8  shows  that  the  prediction  of  laminar  separation, 
reattachment,  transition  and  separation  points  on  the  airfoil 
surfaces  are  in  reasonable  agreement  with  the  experimental 
data  taken  from  Reference  9. 
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5.8.  Transition  and  Separation  Positions  for  NACA 


VI.  SUMMARY  AND  RECOMMENDATIONS 


I 
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This  thesis  treats  the  problem  of  incompressible  two- 
dimensional  steady  flow  past  airfoils  or  airfoil  combinations 
at  large  angles  of  attack.  A  panel  method  was  developed  to 
compute  the  inviscid  flow  over  two  cylinders,  airfoil-flap 
combinations  and  airfoils  in  ground  effect.  In  addition, 
Cebeci's  viscous/inviscid  interaction  method  was  applied  to 
several  airfoils'  and  compared  |^v|th  available  experimental 
data.  The  agreement  is  generally  quite  encouraging.  The 
calculations  show  the  sensitivity  of  the  computations  to 
Reynolds  number  and  transition.  More  work  is  required  to 
evaluate  the  potential  and  limitations  of  the  viscous/inviscid 
interaction  method. 
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